#--------------------------------------------------#
#   Academic freedom and the onset of autocratization   #
#--------------------------------------------------#

# Title: "Academic freedom and the onset of autocratization" #
# Authors: "Pelke, Lars", FAU Erlangen-Nürnberg
# date: 2023-04-24
# journal: Democratization
# DOI: 10.1080/13510347.2023.2207213
# written under "R version 4.2.3 (2023-03-15)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# set working directory
#setwd("P:/PIPM/7. VWS Projekt Academic Freedom Index/replication_files_autocratization_and_academic_freedom")

#### Libraries ####

library(tidyverse)
library(performance)
library(estimatr)
library(ggeffects)
library(ggokabeito)

library(marginaleffects)
library(ggpubr)

#### Load data ####

merged_data <- readRDS("data/merged_data_master.rds") 

merged_data_full <- merged_data %>%
  drop_na(v2x_polyarchy, v2xca_academ, e_pop, e_gdppc, v2xca_academ_10)

merged_data_full <- merged_data_full %>%
  drop_na(democratic_values)


#### Delete Cohorts with less than 1000 persons 

merged_data_full %>%
  group_by(cohort_10.x) %>%
  count()

merged_data_full <- merged_data_full %>%
  filter(cohort_10.x >= 1915) # delete all respondents which birth cohort is before 1910

merged_data_full <- merged_data_full %>%
  dplyr::select(-c(country_name.y, cohort_10.y)) %>%
  rename(country_name = country_name.x, 
         cohort_10 = cohort_10.x)

#### Generate country_cohort variable for calculating std. err. ####

merged_data_full <- merged_data_full %>%
  mutate(country_cohort = str_c(country_text_id, cohortmatch10_20, sep = "_"))

table(merged_data_full$country_cohort)

#### Poland Case Study ####

poland_data <- merged_data_full %>%
  filter(country_name == "Poland") %>% 
  group_by(year) %>%
  summarize(mean_support = mean(democratic_values_normalized, na.rm = T), 
            sd_support =  sd(democratic_values_normalized, na.rm = T))

#################################
#Analysis
##################################

## Model 1 without macro-level vars##
m1 <- lm_robust(democratic_values_normalized ~  graduate + sex + age + income_deciles + children +
                  unemployed, 
                fixed_effects = ~ year + country_id + cohortmatch10_20,
                clusters = country_cohort, se_type = "stata", 
                data = merged_data_full)
summary(m1)

## Model 2 ##
m2 <- lm_robust(democratic_values_normalized ~  graduate + sex + age + income_deciles + children +
                  unemployed + v2xca_academ_10 + v2x_polyarchy + log10(e_gdppc) + log10(e_pop) , 
                fixed_effects = ~ year + country_id + cohortmatch10_20,
                clusters = country_cohort, se_type = "stata", 
                data = merged_data_full)
summary(m2)


## Model 3 ##
m3 <- lm_robust(democratic_values_normalized ~  graduate + sex + age + income_deciles + children +
                  unemployed + v2xca_academ_10 + v2x_polyarchy + log10(e_gdppc) + log10(e_pop) + 
                  graduate:v2xca_academ_10 + as.factor(year) + as.factor(cohortmatch10_20), 
                fixed_effects = ~ country_id,
                clusters = country_cohort, se_type = "stata", 
                data = merged_data_full)
summary(m3)

options(modelsummary_factory_html = "kableExtra")

modelsummary::modelsummary(list(m1, m2, m3),
                           output = "results/Table1_10.tex", 
                           fmt = 3,
                           estimate  = c("{estimate}{stars}"),
                           statistic = 'conf.int', 
                           conf_level = .95, 
                           coef_omit = "year|country_id|cohortmatch10_20", 
                           coef_rename = c("graduate" = "Graduate", 
                                           "sex" = "Sex", 
                                           "age" = "Age", 
                                           "income_deciles" = "Income (ten cat.)", 
                                           "children" = "Children", 
                                           "unemployed" = "Unemployed", 
                                           "v2xca_academ_10" = "Academic_Freedom at C", 
                                           "v2x_polyarchy" = "Electoral Democracy Score", 
                                           "log10(e_gdppc)" = "GDP pc(log)", 
                                           "log10(e_pop)" = "Population (log)", 
                                           "graduate:v2xca_academ_10" = "Academic Freedom at C * Graduate"))


modelsummary::modelsummary(list(m1, m2, m3),
                           output = "results/Table1_10.docx", 
                           fmt = 3,
                           estimate  = c("{estimate}{stars}"),
                           statistic = 'conf.int', 
                           conf_level = .95, 
                           coef_omit = "year|country_id|cohortmatch10_20", 
                           coef_rename = c("graduate" = "Graduate", 
                                           "sex" = "Sex", 
                                           "age" = "Age", 
                                           "income_deciles" = "Income (ten cat.)", 
                                           "children" = "Children", 
                                           "unemployed" = "Unemployed", 
                                           "v2xca_academ_10" = "Academic_Freedom at C", 
                                           "v2x_polyarchy" = "Electoral Democracy Score", 
                                           "log10(e_gdppc)" = "GDP pc(log)", 
                                           "log10(e_pop)" = "Population (log)", 
                                           "graduate:v2xca_academ_10" = "Academic Freedom at C * Graduate"))


#### Figure 2: Coefficient Plot ####


b <- list(geom_vline(xintercept = 0, color = 'orange'),
          annotate("rect", alpha = .1,
                   xmin = -.5, xmax = .5, 
                   ymin = -Inf, ymax = Inf),
          geom_point(aes(y = term, x = estimate), alpha = .3, 
                     size = 10, color = 'red'))

modelsummary::modelplot(m3, coef_omit = 'Intercept|year|country_id|cohortmatch10_20', background = b)



##Simulations for data 

merged_data_full$cohortmatch10_20 <- as.factor(merged_data_full$cohortmatch10_20)
merged_data_full$year <- as.factor(merged_data_full$year)



## Model 3 ##
m3_sim <- lm(democratic_values_normalized ~  graduate + sex + age + income_deciles + children +
               unemployed + v2xca_academ_10 + v2x_polyarchy + log10(e_gdppc) + log10(e_pop) + 
               graduate:v2xca_academ_10 + year + cohortmatch10_20 + country_text_id, 
             data = merged_data_full)
summary(m3_sim)


m3_pred <- ggpredict(m3_sim, terms = c("v2xca_academ_10[all]", "graduate[0,1]"),
                     type = "fe", condition = c(cohortmatch10_20 = "1965"), 
                     cov.fun = "vcovCL", 
                     vcov.type = "HC1",
                     vcov.args = list(cluster = merged_data_full$country_cohort))

f1b <- ggplot(m3_pred, aes(x = x, y = predicted, group = group, color = group, fill = group)) +
  geom_line() +
  geom_ribbon(aes(ymin= conf.low, ymax = conf.high), alpha = 0.3) +
  theme_bw() +
  ggtitle("") +
  labs(color= "Graduate", fill = "Graduate",  
       y= "Predicted Democratic Support",
       x = "Academic Freedom at University Socialization") + 
  theme(legend.position="bottom") +
  scale_color_okabe_ito(order = c(2, 6)) +
  scale_fill_okabe_ito(order = c(2, 6)) 


#### Predict Marginal effects ####

margin_m3 <- marginaleffects(m3_sim,  variables = "v2xca_academ_10", vcov = ~country_cohort,
                        newdata = datagrid(graduate = c(0,1)))

margin_m3$graduate <- as.factor(margin_m3$graduate)
levels(margin_m3$graduate) <- c("Non-Graduate", "Graduate")

f1a <- ggplot(margin_m3, aes(x = graduate, y = dydx, group = graduate, color = graduate, fill = graduate)) +
  geom_point() +
  geom_linerange(aes(ymin = conf.low, ymax=  conf.high)) + 
  theme_bw() +
  ggtitle("") +
  labs(color= "Graduate", fill = "Graduate",  
       y= "Average Marinal Effect of Academic Freedom at C",
       x = "Graduate") + 
  theme(legend.position="bottom") +
  scale_color_okabe_ito(order = c(2, 6)) +
  scale_fill_okabe_ito(order = c(2, 6)) 


ggarrange(f1a, f1b)


ggsave("figures/Figure_2_10.pdf", width = 25,
       height = 15,
       units = c("cm"), dpi = 1000)


###########################################

m3_pred <- ggpredict(m3_sim, terms = c("v2xca_academ_10[all]", "graduate[0,1]", "cohortmatch10_20[all]"),
                     type = "fe", 
                     cov.fun = "vcovCL", 
                     vcov.type = "HC1",
                     vcov.args = list(cluster = merged_data_full$country_cohort))

f2 <- ggplot(m3_pred, aes(x = x, y = predicted, group = group, color = group, fill = group)) +
  geom_line() +
  geom_ribbon(aes(ymin= conf.low, ymax = conf.high), alpha = 0.3) +
  theme_bw() +
  ggtitle("") +
  facet_wrap(vars(facet), nrow = 4) +
  labs(color= "Graduate", fill = "Graduate",  
       y= "Predicted Democratic Support",
       x = "Academic Freedom at University Socialization") + 
  theme(legend.position="bottom") +
  scale_color_okabe_ito(order = c(2, 6)) +
  scale_fill_okabe_ito(order = c(2, 6)) 


ggsave(plot = f2, "figures/Figure_3_10.pdf", width = 20,
       height = 20,
       units = c("cm"), dpi = 1000)





